\section{Desarrollo}

En esta sección detallaremos los métodos implementados para resolver los puntos pedidos en el trabajo práctico.

Antes de comenzar cabe aclarar que a pesar de que las imágenes suministradas son imágenes en formato RGB, como puede suponerse que 
la superficie de los objetos fotografiados es lambertiana, estas fueron convertidas a escala de grises. Es por 
ello que al referirnos al \quotel valor\quoter de un píxel, nos referimos al valor numérico del píxel en escala de grises, 
el cual tiene un rango numérico de 0 a 255.

\subsection{Calibración}

Como primer punto para la resolución, necesitamos obtener las direcciones de iluminación de cada imágen. Esto es, la dirección 
desde la cual está iluminada cada imágen a procesar. Esta dirección puede representarse como un vector de 3 dimensiones.

Como podemos suponer que la fuente de iluminación está lo suficientemente alejada del objeto, el vector de iluminación será 
único para todos los puntos de la imágen, es por ello que al conseguir la dirección de iluminación en el punto más brillante 
de la imágen, esa dirección puede utilizarse para el resto de los puntos.

\subsubsection{Parseo de la Esfera}

Como podemos suponer que las imágenes que tienen el mísmo número están iluminadas del mismo ángulo, para la obtención de la 
dirección de iluminación utilizamos las imágenes de las esferas respectivas a las imágenes reales a procesar. Esto se 
debe principalmente a que en la esfera la dirección de iluminación será normal a su superficie, y las direcciones normales 
a una esfera son fáciles de calcular. Es por esto que necesitamos parsear las imágenes de la esfera.

Para hacerlo, simplemente levantamos los valores de los píxeles de cada imágen y los guardamos dentro de una matriz.

\subsubsection{Centrar la esfera y obtener radio}

Para que las direcciones de iluminación tengan sentido, el eje de coordenadas de las mismas debe estar ubicado en el centro geométrico de 
la esfera. Es por ello que debemos encontrarlo, para relativizar todas las coordenadas dentro de las imágenes al centro 
de la misma.

Para hacerlo utilizamos el archivo de la máscara de la esfera. Hacer esto es sencillo: encontramos el punto perteneciente 
a la esfera ubicado más arriba en la imágen (llamémoslo $sup = (sup_x, sup_y)$) y el ubicado más a la izquierda 
(llamémoslo $izq = (izq_x, izq_y)$). Por ser una esfera, el centro de la misma se encuentra en la fila $izq_y$ y en la columna $sup_x$.

Una vez obtenido el centro de la esfera, encontrar el rádio de la misma también es sencillo. Simplemente contamos la cantidad de 
píxeles que hay entre el centro de la esfera y, por ejemplo, el punto más a la izquierda de la misma.

\subsubsection{Obtención del punto de mayor iluminación} \label{obtencion_iluminacion}

En un primer momento creimos que para obtener el punto de mayor iluminación simplemente debiamos recorrer la imágen de la esfera 
y quedarnos con la coordenada del punto de mayor valor. Pero esta idea no tiene en cuenta el hecho de que pueden existir más de 
dos puntos con el mismo valor máximo ni que pueden existir errores en la medición de las intensidades de los píxeles. Es por 
ello que decidimos descartar esta idea y utilizar una más sofisticada.

La imágen es recorrida en busca de todos los puntos de mayor valor. Una vez encontrados estos $n$ puntos, analizamos los píxeles 
que los rodean. Este análisis consiste en desplazarnos, por cada uno de los $n$ píxeles, en las 4 direcciones en busca de 4 
posiciones $izq$, $der$, $sup$ e $inf$. Nos desplazamos en cada dirección mientras el valor del píxel en el que nos encontramos 
esté un 10\% por debajo del valor máximo. Podemos ver en las figuras \ref{pixel1} y \ref{pixel2} una representación de esta idea:

\begin{figure}[h]
\hfill
\begin{minipage}[t]{.45\textwidth}
\begin{center}
	\includegraphics[scale = 0.30]{graficos/pixel1.jpg}
	\caption{Píxel de mayor valor.}
	\label{pixel1}
\end{center}
\end{minipage}
\hfill
\begin{minipage}[t]{.45\textwidth}
\begin{center}
	\includegraphics[scale = 0.30]{graficos/pixel2.jpg}
	\caption{Posibles píxeles $izq$, $der$, $sup$ e $inf$.}
	\label{pixel2}
\end{center}
\end{minipage}
\hfill
\end{figure}

Una vez obtenidos esas 4 posiciones para cada uno de los $n$ puntos, obtenemos el punto medio $pm_x$ entre las posiciones $izq$ y $der$, 
y el punto medio $pm_y$ entre las posiciones $sup$ e $inf$, como puede verse en la figura \ref{pixel3}. 
Estos dos puntos medios calculados determinan un nuevo punto ubicado en $(pm_x, pm_y)$, como lo indica la figura \ref{pixel4}.

\begin{figure}[h]
\hfill
\begin{minipage}[t]{.45\textwidth}
\begin{center}
	\includegraphics[scale = 0.30]{graficos/pixel3.jpg}
	\caption{Pixeles $pm_x$ y $pm_y$.}
	\label{pixel3}
\end{center}
\end{minipage}
\hfill
\begin{minipage}[t]{.45\textwidth}
\begin{center}
	\includegraphics[scale = 0.30]{graficos/pixel4.jpg}
	\caption{Píxel de mayor intensidad.}
	\label{pixel4}
\end{center}
\end{minipage}
\hfill
\end{figure}

De esta forma cada uno de los $n$ puntos determina un nuevo punto. Sean $n'$ estos puntos nuevos.

Elegimos al punto de mayor intensidad de la imágen como el punto de mayor valor entre los $n'$ puntos. En caso de existir dos o mas puntos con 
el mismo valor, nos quedamos con el primero de ellos. Denominamos $max = (max_x, max_y)$ al punto de mayor iluminación de cada imágen.

\subsubsection{Obtención del vector de iluminación}

Una vez determinada la posición $max$ de mayor iluminación esta debe expresarse en los nuevos ejes de coordenadas encontrados. 
Sea $(cen_x, cen_y)$ la posición del centro de la esfera, la posición $(max_x - cen_x, max_y - cen_y) = (dir_x, dir_y)$ corresponde a la posición del 
píxel de mayor intensidad en el nuevo eje coordenadas. Sólo nos resta conocer la \quotel profundidad\quoter que hay en la esfera en dicho punto 
para obtener el vector de dirección de iluminación. Conociendo el rádio de la esfera, nuevamente es sencillo. Siendo $r$ el radio de 
la esfera y $dir_z$ la profundidad buscada, sabemos que:

\begin{center}
$r^2 = {dir_x}^2 + {dir_y}^2 + {dir_z}^2$
\end{center}

Luego podemos obtener la profundidad de la siguiente forma:

\begin{center}
$dir_z = \sqrt{r^2 - {dir_x}^2 - {dir_y}^2}$
\end{center}

De esta manera conocemos las tres componentes que conforman al vector de iluminación de una imágen en particular: $v = (dir_x, dir_y, dir_z)$.

\subsection{Campo Normal}

Una vez conseguidas las direcciones de iluminación de las imágenes, el siguiente paso es construir el campo normal de la figura en 
cuestión.

\subsubsection{Parseo y recorte de imágenes} \label{parrec}

En una primera aproximación realizamos esto procesando la totalidad de los píxeles de las imágenes. Esto es una mala idea, ya que 
en comparación gran cantidad de píxeles pertenecen al fondo de las mismas, y se pardía tiempo y espacio en almacenar y procesar 
algo que no era necesario hacerlo. Es por ello que decidimos recortar las imágenes utilizando sus máscaras. Este recorte 
consiste en quedarnos sólo con en menor rectángulo en el cual quepa la imágen.

Para realizar esto parseamos la máscara del objeto a procesar, guardándola en una matriz y determinando luego 
las dos filas $f_1$ (por arriba) y $f_2$ (por abajo) y dos columnas $c_1$ (por la izquierda) y $c_2$ (por la derecha) 
que contienen al objeto de la imágen con mayor exactitud.

Una vez hecho esto parseamos las tres imágenes a procesar. Estas imágenes son almadenadas dentro de una matriz, con $f_2$ - $f_1$ filas 
y $c_2$ - $c_1$ columnas. A pesar de que las imágenes son recortadas, una parte del fondo es almacenada dentro de la matriz. Los 
valores de estos píxeles en las imágenes originales no son exactamente cero, aunque si un valor muy cercano. Ya que conceptualmente estos 
valores son cero, cuando almacenamos en la matriz un píxel perteneciente al fondo de la imágen, guardamos un cero. En cambio, los valores 
pertenecientes al objeto son guardados con el mismo valor que tiene en la imágen.

\subsubsection{Construcción del campo normal}

Para obtener el campo normal de la figura necesitamos resolver un sistema de ecuaciones lineales de tres ecuaciones con tres incógnitas, que 
corresponde con la ecuación (5) en el enunciado del trabajo práctico, el cual puede encontrarse en la sección \ref{enunciado}.

Esta resolución es llevada a cabo por la clase \textit{sistema} (implementada en el archivo \texttt{sistema.cpp}). 
En ella se triangula la matriz del sistema aplicando el método de Gauss con pivoteo parcial 
y despejando las incógnitas utilizando \textit{backward substitution}.

Disponiendo entonces de una forma para resolver sistemas de ecuaciones lineales, procedemos a construir el campo normal de la figura. 
Para ello recorremos los píxeles de cada una 
de las tres imágenes que tenemos almacenadas en matrices. Para cada uno de ellos, armamos el siguiente sistema de ecuaciones:

\begin{equation*}
\begin{bmatrix}
I_1\\
I_2\\
I_3
\end{bmatrix}
=
\begin{bmatrix}
s_{x}^{1} & s_{y}^{1} & s_{z}^{1}\\
s_{x}^{2} & s_{y}^{2} & s_{z}^{2}\\
s_{x}^{3} & s_{y}^{3} & s_{z}^{3}
\end{bmatrix}
\begin{bmatrix}
m_x \\
m_y \\
m_z 
\end{bmatrix}
\end{equation*}

Donde $I_i$ representa la intensidad del píxel que está procesándose en la i-ésima imágen, $s_{x}^{i}$, $s_{y}^{i}$ y $s_{z}^{i}$ 
representan las componentes $x$, $y$ y $z$ 
del vector de iluminación para la i-ésima imágen. Este sistema es resuelto para cada uno de los píxeles de la imágen. 

Una vez obtenida la solución $(m_x, m_y, m_z)$, esta es normalizada y almacenada en una matriz de triplas, con las mismas dimensiones 
que las matrices de las imágenes recortadas. Sea $m$ dicha matriz, y sea $img^i$ la i-ésima matriz que contiene los valores de los píxeles 
de la i-ésima imágen recortadas, $m_{f, c}$ corresponderá al vector normal a la superficie del objeto para el píxel $img^{i}_{f, c}$.

Es posible que la norma de este vector sea cero. En tal caso no puede ser normalizado. 
Esto sucede en general con los píxeles correspondientes al fondo de la imágen, ya que los valores de $I_1$, $I_2$ y $I_3$ 
son iguales a cero (pues fueron seteados en ese valor como se explica en la sección \ref{parrec}). 
De encontrarnos con un vector con norma cero, el vector almacenado en $m_{f, c}$ es $(0, 0, -1)$. Fue elegido este valor 
ya que tiene norma 1 y podrá ser diferenciado del resto de los vectores normales.

De no tratarse de un vector con norma cero, el vector es normalizado y almacenado en $m_{f, c}$. Una vez procesados todos los píxeles 
de la imágen de esta forma, obtenemos el campo normal en la matriz $m$.

\subsection{Hallar la profundidad}

Una vez determinado el campo normal de la imágen, el siguiente paso es determinar la profundidad de cada píxel.

\subsubsection{Obtención del sistema $M x = v$}

Siendo $f$ y $c$ la cantidad de filas y columnas respectivamente de la matriz $m$, nos basamos en las ecuaciones 
(11) y (12) del enunciado (Sección \ref{enunciado}) para construir un sistema de ecuaciones de la siguiente forma:

\begin{center}
$M x = v$ 
\end{center}

Siendo $M$ una matriz de $(c * f * 2)$ filas y $((f + 1) * (c + 1) - 1)$ columnas, $x$ un vector columna 
de $((f + 1) * (c + 1) - 1)$ filas y $v$ un vector columna de $(c * f * 2)$ filas. Explicaremos a 
continuación con más detalle la forma de generar la matriz $M$ y el vector $v$.

Ya que las ecuaciones (11) y (12) utilizan los valores de los píxeles ubicados a la derecha y debajo 
de un determinado píxel, estas ecuaciones no están bien definidas para los píxeles de la última fila 
y última columna de la matriz $m$. Es por ello que debemos pensar que la matriz $m$ contiene una fila 
y una columna adicional, como lo muestra la Figura \ref{matrizext}.

De esta forma, la cantidad de incógnitas a despejar es igual a $((f + 1) * (c + 1) - 1)$. 
Es por ello que la cantidad de columnas de $M$ es también $((f + 1) * (c + 1) - 1)$. 
Además, como por cada píxel se tienen dos ecuaciones y hay $f * c$ píxeles, la cantidad 
de filas de $M$ es $(2 * f * c)$.

Nuevamente, en base a las ecuaciones (11) y (12), siendo $m_{i, j} = (x_{i, j}\ ,\ y_{i, j}\ ,\ z_{i, j})$, 
armamos el término independiente $v$ de la siguiente forma:

\begin{center}
$v^t = (-x_{1, 1}\ ,\ -y_{1, 1}\ ,\ -x_{1, 2}\ ,\ -y_{1, 2}\ ,\ \hdots\ ,\ -x_{1, c}\ ,\ -y_{1, c}\ ,\ -x_{2, 1}\ ,\ -y_{2, 1}\ ,\ \hdots\ ,\ -x_{f, c}\ ,\ -y_{f, c})$
\end{center}

Notar que los píxeles de la imágen son recorridos por filas, agregando al vector primero 
la componente $x$, luego la $y$, ambas multiplicadas por $(-1)$, para luego pasar al siguiente píxel y repetir el proceso.

\begin{figure}[h] % Acomodar este donde mejor quede
\hfill
\begin{center}
	\includegraphics[scale = 0.7]{graficos/imgext.png}
	\caption{Matriz de píxeles extendida.}
	\label{matrizext}
\end{center}
\hfill
\end{figure}

La forma en la que decidimos armar el vector $v$ determinó la forma que debía tener la matriz $M$ para que formen las 
ecuaciones (11) y (12) del enunciado para cada píxel: La primera fila corresponde a los coeficientes de la ecuación 
(11) para el píxel $m_{1,1}$, la segunda fila corresponde a los coeficientes de la ecuación (12) para el píxel 
$m_{1,1}$, la tercera fila corresponde a los coeficientes de la ecuación (11) para el píxel $m_{1,2}$, la cuarta 
fila corresponde a los coeficientes de la ecuación (12) para el píxel $m_{1,2}$, etcétera. Notar que al armar la matriz 
de esta forma, esta contiene exactamente dos elementos distintos de cero por fila: La gran mayoría de las posiciones 
de la matriz contienen ceros, es decir, es una matriz esparsa.

Las columnas de $M$ corresponden a las incógnitas buscadas, es decir, las profundidades de los píxeles. Sea $z_{i,j}$ 
la profundidad del píxel $m_{i,j}$, las columnas de $M$ están ordenadas de forma tal que el vector de incógnitas 
$x$ es:

\begin{center}
$x^t = (z_{1,1}\ ,\ z_{1,2}\ ,\ \hdots\ ,\ z_{1,c}\ ,\ z_{1,c+1}\ ,\ z_{2,1}\ ,\ \hdots\ ,\ z_{f, 1}\ ,\ \hdots\ ,\ z_{f,c}\ ,\ 
z_{f, c+1}\ ,\ z_{f+1, 1}\ ,\ \hdots\ ,\ z_{f+1, c})$
\end{center}

Notar que este vector contiene más incógnitas que las buscadas. Esto se debe a los píxeles adicionales que fueron 
\quotel agregados\quoter 
a la matriz $m$. Estas incógnitas \quotel extras\quoter son aquellas que tienen como columna el valor $c+1$ y/o como fila el valor $f+1$. 
Estas deberán ignorarse al momento de recostruir la representación en 3D del objeto.

%Mostramos a continuación un pequeño ejemplo de una matriz $M$ contruida en base a una matriz $m$ de 2 * 2.

%\[\begin{bmatrix}
%-n_{11}&\vec{0} & n_{11}&0&0&0&0\\
%-n_{11}&n_{11}&0&0&0&0&0&0\\
%0&\ -n_{12} &\vec{0}& n_{11}&0&0&0 \\
%0&-n_{12}&n_{12}&0&0&0&0&0\\
%0&0&0&-n_{21} &\vec{0} & n_{21}&0 \\
%0&0&0&-n_{21}&n_{21}&0&0&0 \\
%0&0&0&0&-n_{22}&\vec{0} & n_{22}\\
%0&0&0&0&-n_{22}&n_{22}&0&0\\ 

%\end{bmatrix}\]
%En este caso $\vec{0}$ es un vector (0,0), pero en el caso mas general de la imágen es en realidad de la dimension de
%la cantidad de columnas de la imágen.

%El patrón nos permite ver que cada 2 veces la cantidad de filas de la imágen original hay un desplazamiento extra
%en la cantidad de ceros y luego el patrón se repite.

%\[\begin{bmatrix}
%-n_{11}&\vec{0} & n_{11}&0&\cdots &0 \\
%-n_{11}&n_{11}&0&\cdots&0 \\
%0&\ -n_{12} &\vec{0}& n_{11}&0& \cdots &0 \\
%0&-n_{12}&n_{12}&0&\cdots&0 \\
%\vdots & \vdots & \ddots & \vdots\\
%\vec{Vcantcols-1}&-n_{1cols+1}&\vec{0} & n_{1cols+1}&0&\cdots &0 \\
%0&\ n_{11} &&\cdots &0 \\
%0&-n_{12}&n_{12}&0&\cdots&0 \\

%a_{n1}&a_{n2}&\cdots &a_{nn}
%\end{bmatrix}\]


\subsubsection{Implementación de la matriz esparsa}

Debido al gran tamaño de la matriz $M$ una implementación clásica de la misma 
(como por ejemplo un vectores de vectores), no es factible. No solo por la 
cantidad de memoria necesaria para almacenarla, sino por el tiempo que tomaría 
aplicarle distintos algoritmos, como por ejemplo de triangulación. Es por ello que 
esta idea fue descartada de inmediato.

Pero este problema puede solucionarse aprovechando el hecho de que la matriz es 
esparsa. De esta forma, puede tenerse una estructura que represente la matriz, pero 
que solo almacene a las posiciones que contienen elementos distintos de cero, que 
en relación con la cantidad de elementos son muy pocos.

Nuestra primera idea fue utilizar un vector de vectores de tuplas: Un vector principal 
del tamaño de la cantidad de filas totales de la matriz, el cual contiene tuplas que 
representan a las columnas de $M$. Estas tuplas contienen un entero que representa la 
columna a la cual pertenece el elemento y el elemento en si (un double), el cual debe ser distinto de cero. 
Esta estructura se corresponde a lo que denominamos una representación \quotel por filas\quoter de la matriz. 
Esta representación proporciona un rápido acceso a las filas de la matriz, pero encontrar un 
elemento en una columna particular es lineal o logarítmico, ya que los elementos en el vector 
están ordenados en orden ascendente por columnas. Puede verse en la Figura \ref{matrizesp} una 
representación gráfica de esta estructura.

\begin{figure}[h] % Acomodar este donde mejor quede
\hfill
\begin{center}
	\includegraphics[scale = 0.8]{graficos/matriz_esp.png}
	\caption{Matriz esparsa, representación por filas.}
	\label{matrizesp}
\end{center}
\hfill
\end{figure}

Análoga a esta representación por filas, tenemos la representación \quotel por columnas\quoter. 
La estructura no cambia, salvo por la forma de interpretarla. En este caso, el vector 
principal representa las columnas de la matriz y el primer elemento de cada nodo representa la fila.

Para determinar el valor de una posición de la matriz basta con acceder a la fila (columna) 
correspondiente y luego buscar el elemento que tenga la columna (fila) buscada. En caso de no 
encontrarse, significa que el elemento vale cero.

Notar que es deseable que, de trenerse estas dos representaciones para una misma matriz, estas 
sean consistentes entre si. En otras palabras, que ambas representen a la misma matriz.

El archivo \texttt{matriz\_esp.cpp} contiene la implementación de esta matriz 
junto con sus métodos. De todas formas en las secciones siguientes explicaremos con mayor 
detalle algoritmos tales como la multiplicación de $M^t * M$, su triangulación y el despeje de incógnitas.

\subsubsection{Multiplicación $M^t$ por $M$}

Como ya hemos mencionado el sistema a resolver es el siguiente:

\begin{center}
$M x = v$
\end{center}

Ya que la cantidad de filas de $M$ es mucho mayor que la cantidad de columnas, este sistema está sobredeterminado. 
Por esto decidimos utilizar la ecuacion (14) del enunciado y generar el siguiente sistema aproximado sobre el cual operar:

\begin{center}
$M^t * M * x = M^t * v$ 
\end{center}

De esta forma el sistema a resolver sería cuadrado, los cuales son mucho más sencillos algoritmicamente. 
El método que empleamos para su resolución es triangularlo mediante el algoritmo de \textit{Gauss-Jordan} 
y luego despejar las incógnitas utilizando \textit{backward substitution}.

Para generar la matriz $M^t * M$ (que también será esparsa) utilizamos el siguiente algoritmo:

\hfill
\begin{verbatim}
            para entero i = 0, mientras i < filas(M), hacer
            
                  par <entero, double> par_elem_1 = fila(M, i)[0]
                  par <entero, double> par_elem_2 = fila(M, i)[1]
    
                  entero elem_1_colum = first(par_elem_1)
                  entero elem_2_colum = first(par_elem_2)
    
                  double elem_1 = second(par_elem_1)
                  double elem_2 = second(par_elem_2)

                  mtxm[elem_1_colum][elem_1_colum] += elem_1 * elem_1
                  mtxm[elem_1_colum][elem_2_colum] += elem_1 * elem_2
                  mtxm[elem_2_colum][elem_1_colum] += elem_2 * elem_1
                  mtxm[elem_2_colum][elem_2_colum] += elem_2 * elem_2
    
                  i = i + 1          
            
            fin para        
\end{verbatim}
\hfill

Notar que multiplicar la traspuesta de una matriz $m$ por dicha matriz resulta en la multiplicación de las columnas de $m$ entre sí.

Decidimos aprovechar el hecho de que $M$ tiene dos elementos no nulos por fila. Por ello el algoritmo 
utizado no es el convencional para la multiplicación de matrices, sino que es una suerte de \quotel 
multiplicaciones parciales\quoter.

Utilizamos la representación por filas de $M$ para recorrer las filas una a una: Sean $a$ y $b$ 
los únicos dos elementos de una fila de $M$ y sean $c_a$ y $c_b$ las columnas donde se encuentran 
estos elementos, en algún momento se deberán multiplicar $c_a$ con $c_a$, $c_a$ con $c_b$, $c_b$ 
con $c_a$ y $c_b$ con $c_b$. Es por ello que, como $a$ y $b$ están en la misma fila, se deberán 
hacer estos cuatro productos anteriores correspondientes entre $a$ y $b$ almacenando el 
resultado donde corresponda, dentro de una nueva matriz esparsa.

Todas las filas de $M$ son recorridas aplicando estos productos y al finalizar tenemos 
la matriz esparsa que representa a $M^tM$.

Por lo tanto pudimos resolver la multiplicación de matrices en tiempo lineal con respecto 
a la cantidad de elementos de la matriz, que en este caso es igual en ambas matrices 
porque trabajamos con M y M$^t$ que tienen igual dimensión y los mismos elementos.

\subsubsection{Multiplicación de $M^t$ por $v$}

La siguiente tarea a realizar es la multiplicación de $M^t$ por el vector $v$. 
Para hacer esto utilizamos una idea muy similar a la multiplicación de $M^t * M$.

Sean $a$ y $b$ los dos elementos de la fila número $i$ de $M$ y sean $c_a$ y $c_b$ las columnas 
sobre la que están $a$ y $b$, respectivamente, y sea $v' = M^t * v$. 
Multiplicamos el elemento $a$ por $v_i$ y acumulamos el resultado en $v'_{[c_a]}$ y 
multiplicamos al elemento $b$ por $v_i$ y acumulamos el resultado en $v'_{[c_b]}$.     

Al recorrer todas las filas de $M$ realizando este procedimiento obtenemos $M^t * v$.

\subsubsection{Triangulación del sistema $M^t M * x = M^t v$}

Como ya fue mencionado, para resolver este sistema utilizaremos el método de Gauss-Jordan para triangular la matriz $M^t M$. 
Nuevamente, no puede aplicarlse el algoritmo de forma ingenua. Es necesario adaptarlo a la estructura de la matriz para 
aprovechar sus propiedades.

Para la triangulación utilizaremos la representación por columnas de $M^t M$ (a partir de ahora $\tilde{M}$) ya que necesetaremos anular elementos en cada 
columna y de esta forma podemos acceder a ellos rápidamente. Detallamos a continuación el algoritmo:

Durante la explicación se hará referencia a $\tilde{M}_{ij}$, que corresponde al elemento ubicado en la 
fila $i$ columna $j$ de $\tilde{M}$. También se hará referencia a $\tilde{M}_i$, que corresponde a la 
$i$-ésima columna de $\tilde{M}$ (recordar que utilizaremos la representación por columnas).

Comenzamos un ciclo donde recorremos las columnas de $M^t M$. 
Sea $c$ la columna actual a triangular de $\tilde{M}$ realizamos los siguientes pasos:

\begin{itemize}\itemsep2pt

\item Obtenemos el elemento $\tilde{M}_{c,c}$ de la diagonal. Para hacer esto recorremos el vector $\tilde{M}_{c}$ desde 
el comienzo al final viendo la fila en la que se encuentra cada elemento.

Sea $f$ el valor de la fila del elemento que estamos viendo, si $f$ es menor que $c$, pasamos al siguiente elemento ya que 
todavía no llegamos a la diagonal. Si $f$ es igual a $c$ tenemos el elemento de la diagonal.

No podría suceder que no encontremos el elemento de la diagonal. Si esto ocurriese el elemento sería cero. 
Esto se debe a que las filas que COMPLETAR.

\item Ahora calculamos los valores de los multiplicadores para cada fila. Recorremos los elementos de la columna $c$ 
desde el final hacia el principio hasta llegar a la diagonal, guardando en un vector el valor del multiplicador 
junto con la fila a la que está asociado. El vector es recorrido en sentido inverso ya que sabemos que todos los 
elementos que están por debajo de la diagonal deben ser anulados, pero no sabemos en qué posición del vector está 
el elemento diagonal.

Ya que los elementos que están por debajo de la diagonal en la columna $c$ queremos que sean cero, luego de calcular el 
multiplicador asociado a la fila del elemento, este \quotel nodo\quoter es borrado de la columna. Esto equivale 
a ponerlo en 0.

\item Una vez calculados todos los multiplicadores, los aplicamos a las filas correspondientes. 
En vez de recorrer los elementos de las filas en orden, como seria natural, recorremos los elementos de las columnas. 
Recorremos todas las columnas mayores que $c$.

	\begin{itemize}\itemsep4pt

	\item Sea $c'$ la columna actual que estamos recorriendo, debemos encontrar el elemento $\tilde{M}_{cc'}$. Es decir, el elemento 
	que se encuentra en la misma fila que el elemento de la diagonal pero en la columna que estamos procesando. Este es el elemento 
	que acompañará a los multiplicadores en la columna $c'$.

	Para encontrarlo recorremos los elementos de $M_{c'}$ en orden viendo en que filas se encuentran. Sea $f$ la fila en la que se encuentra 
	el elemento actual, si $f$ es menor que $c$ pasamos al siguiente elemento. Si $f$ es igual a $c$, entonces dicho valor es $\tilde{M}_{cc'}$. 
	Si en cambio $f$ es mayor que $c$, como los elementos estan ordenados en orden ascendente por el número de fila a la que pertenece, significa 
	que el elemento buscado es cero.
	
	\item En el caso de que $\tilde{M}_{cc'}$ sea igual a cero, la columna $c'$ no sufrirá modificaciones, por lo que se pasa a la siguiente columna.
	
	\item Si $\tilde{M}_{cc'}$ es distinto de cero, entonces la columna $c'$ debe modificarse. En vez de recorrer la columna $c'$ recorremos 
	el los elementos del vector de multiplicadores, ya que si recorriésemos los elementos de la columna tendríamos que determinar si la fila 
	en la que se encuentra cada elemento tiene asociado un multiplicador (si no lo tiene el elemento no se modifica, pero si lo tiene, sí). 
	En cambio si recorremos los multiplicadores, nos aseguramos que el elemento deberá modificarse sin importar su valor.
	
	Siendo $mult$ el multiplicador actual que estamos viendo, y siendo $f'$ la fila asociada a dicho multiplicador buscamos en 
	$\tilde{M}_{c'}$ el elemento ubicado en la fila $f'$. Una vez encontrado, calculamos el nuevo valor $val = \tilde{M}_{f'c'} + \tilde{M}_{cc'} * mult$ 
	que irá en la posición $\tilde{M}_{f'c'}$.
	
		\begin{itemize}\itemsep8pt
		\item En el caso de que $\tilde{M}_{f'c'}$ sea igual a cero, el nodo correspondiente a este elemento no existía en el vector $\tilde{M}_{c'}$. 
		Luego $val$ es agregado al vector junto con la fila a la cual pertenece.
		
		\item Si $\tilde{M}_{f'c'}$ no era igual a cero, hay dos posibilidades. Si $val$ es igual a cero, entonces hay que eliminar el nodo de 
		$\tilde{M}_{c'}$, ya que no se almacenan ceros. Si $val$ no es igual a cero, simplemente reemplazamos el valor del nodo por $val$.
		\end{itemize}

	\end{itemize}
	
\item Luego de recorrer todas las columnas mayores a $c$, debemos aplicar también los multiplicadores al vector $M^t v$. Procedemos de igual forma que con el resto de las columnas, pero teniendo en cuenta de que $M^t v$ es un vector con todas las posiciones. Obtenemos el elemento ubicado en la posición 
$c$ del vector, que será el que acompañe al multiplicador. Luego recorremos el vector de multiplicadores aplicandolos a las posiciones correspondientes. Con la única diferencia de que el número obtenido es almacenado en el vector sin importar cual sea su valor.
\end{itemize}

Al finalizar el algoritmo, tendremos a la matriz $\tilde{M}$ triangulada sin ceros en la diagonal.

Notar que en este momento la representación por filas de $M^t M$ no es consistente con la representación por columnas. Ya que para el despeje de incógnitas sería cómoda la representación por filas, reestablecemos el invariante de la representación por filas. Para lograrlo simplemente limpiamos la estructura de la representación por filas y la completamos con los valores correspondientes.

\subsubsection{Despeje de incógnitas}

Con la matriz triangulada sin ceros en la diagonal sólo resta despejar las incógnitas. 
Aplicamos el algoritmo de \textit{backward substitution}, adaptado a nuestra estructura.

El algoritmo es aplicado de la forma natural, pero utilizando la representación por filas 
para obtener rápidamente los elementos de cada fila distintos de cero. Estos son buscados 
desde el final de la fila hacia el elemento de la diagonal, el cual siempre existe, 
por lo que es distinto de cero.

Al finalizar el despeje tenemos $(f +1) * (c + 1) - 1$ incógnitas despejadas, 
siendo $f$ y $c$ la cantidad de filas y columnas de las imágenes recortadas.

Como ya fue mencionado, no todas estas incógnitas son las que nos interesan, 
sino que debemos desechar las que corresponden a los píxeles agregados a la matriz $m$ de normales.

Por la forma en la que fue construida $M$, sabemos que cada $c$ elementos en 
el vector $x$ tenemos una de las incógnitas \quotel extras\quoter, al igual 
que las últimas $c+1$ incógnitas.

Estas simplemente son ignoradas al momento de devolver las profundidades buscadas.




